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We develop a systematic method to derive all orders of mode couplings in a weakly nonlinear 
approach to the dynamics of the interface between two immiscible viscous fluids in a Hele-Shaw 
cell. The method is completely general. It includes both the channel geometry driven by gravity 
and pressure, and the radial geometry with arbitrary injection and centrifugal driving. We find the 
finite radius of convergence of the mode-coupling expansion. In the channel geometry we carry out 
the calculation up to third-order couplings, which is necessary to account for the time-dependent 
Saffman-Taylor finger solution and the case of zero viscosity contrast. Both in the channel and the 
radial geometries, the explicit results provide relevant analytical information about the role that the 
viscosity contrast and the surface tension play in the dynamics. We finally check the quantitative 
validity of different orders of approximation against a physically relevant, exact time-dependent 
solution. The agreement between the low order approximations and the exact solution is excellent 
S ■ within the radius of convergence, and reasonably good even beyond that. 

PACS: 47.20.Ma, 47.20.Hw, 47.54.+r, 47.20.Ky 

<N 

The morphological instability of fluid interfaces in Hele-Shaw flows |ljj2| has become a paradigm of interfacial pattern 
formation in nonequilibrium systems . As opposed to the most commonly studied "bulk" pattern forming systems 
Bj the inherent difficulties of the free-boundary problems associated to interfacial growth makes the latter even 
more elusive to analytical treatment. As a prototype of interfacial instabilities in diffusion-limited growth problems 
(including for instance dendritic growth, solidification of mixtures, chemical electrodeposition, flame propagation, 
H ' etc.) the Saffman-Taylor problem R is a relatively simple case, both theoretically and experimentally, well suited to 
^ \ gain insight into generic dynamical features in the broad context of nonlinear interface phenomena. 
0^ ■ The intrinsic difficulty of free-boundary problems is manifest in the fact that the interface dynamics is highly 
non-local. Furthermore, the nature of the instability (except in some cases, such as in directional solidification of 
binary alloys) usually produces non-saturated growth, which inevitably results in a highly nonlinear dynamics. In 
bulk instabilities, when the control parameter is near threshold, the traditional weakly nonlinear techniques lead to 
a universal description of patterns in terms of amplitude equations, based on center manifold reduction |^,^]. These 
techniques, however, are not so useful for interfacial problems in which nonlinearities do not saturate the growth. 
This is the case for instance of viscous fingering. In the case of the channel geometry (Saffman-Taylor problem) [QJ^] 
the interface restabilizcs in a nontrivial morphology (the Saffman-Taylor finger) which keeps growing at a finite rate. 
For circular geometries |j 10 the patterns do not reach an equivalent steady state and the interplay of tip-splitting 
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I. INTRODUCTION 



events and screening effects may result in a variety of complicated morphologies. In these problems all the weakly 
nonlinear techniques apply only to a transient in the early nonlinear regime. Nevertheless, compared to the more 
traditional ones |^|, the weakly nonlinear analysis developed in this paper is not restricted to situations near the 
instability threshold, where a separation of scales is exploited. Instead, we expand on the amplitudes of the whole 
! spectrum of modes. 

In this paper we will deal both with channel and radial geometry Hele-Shaw flows. In the traditional Saffman-Taylor 
problem (channel geometry) the pressure- and gravity-driven instabilities can be formally mapped into each other 
in the appropriate reference frames, so there is really no different interface dynamics for the two physical situations. 
The problem, then, contains two independent dimensionless parameters, namely, a dimensionless surface tension B 
and the viscosity contrast or Atwood ratio A |pd| . In the radial geometry, though, there is no such formal mapping. 
We will see that the injection and the centrifugal forcing are not equivalent and three independent parameters must 
be considered. 

The situation most commonly studied in the literature is the high viscosity contrast limit, A = I , where one of 
the two fluids is non-viscous Q (typically air displacing a viscous fluid). The singular perturbation character of the 
surface tension B has received most of the attention as responsible for the subtle mechanism of steady-state selection, 
namely the fact that surface tension "selects" a single finger solution out of a continuum of solutions for B = [[l2]-[l4j . 
More recently, the crucial role of surface tension in the dynamics of fingering patterns has been pointed out. Siegel 
and Tanveer |I5|-|l7| have shown that the exact, nonsingular time-dependent solutions known for the case with B = 
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may differ significantly from the corresponding ones with B — > + after a time which is of order one (B ). In practice, 
this implies that exact solutions of the problem with B = (including those with no finite time singularities) may lead 
to completely incorrect asymptotic behaviour as compared to the regularized ones, with B -c 1. A careful analysis 
of these questions may be found in Ref. fll8| . Notice however that such analysis restricts to small B, while in many 
cases (for instance, for fingers emerging naturally from the linear instability, with the characteristic length scale of 
the linearly most unstable mode) the effective dimensionlcss surface tension is necessarily B ~ 1. Understanding the 
dynamics of finger competition in typical experimental conditions thus requires considering relatively large values of 
B, for which the perturbative techniques of |15|-|l7|| fail. 

On the other hand, an important role of viscosity contrast A in the dynamics of finger competition has been observed 
both numerically and experimentally |2fj|-p2]|. A careful characterization of the interface evolution has shown 

that for A = the finger competition process is ineffective, and that the system does not approach the usual single 
finger Saffman-Taylor attractor |^,|4| . Although the nature or existence of other attractors is still an open question, 
it seems that the basin of attraction of the Saffman-Taylor solution does depend on A, and is particularly sensitive 
to A in the neighborhood of A ~ 1. In any case, it is clear that the viscosity contrast plays also a crucial role in 
the highly nonlinear regime, and that tuning A in its full range is necessary to elucidate some of the important open 
questions. 

Finally, an interesting interplay between B and A in connection with the selection problem is apparent in that, 
despite the fact that single finger stationary solutions of any width do exist for B = regardless of viscosity contrast 
A, the only single finger time dependent solution of the A = 1 (B = 0) problem which is also a solution for any 
viscosity contrast A is the one that fills one half of the channel [G5p9], which is precisely the solution selected by 
surface tension in the limit B — > 0. Whether deeper consequences concerning the selection problem can be drawn 
from this fact is also an interesting open question. 

In order to gain analytical insight into these dynamical questions we propose here a systematic weakly nonlinear 
expansion of the problem of viscous fingering in Hele-Shaw flows, including all traditional setups and the most recent 
one of rotating flows [ p7| . The basic motivation is to be able to extract information which is nonperturbative in any 
of the two basic parameters, which are taken as completely arbitrary. The expansion parameter will be basically the 
mode amplitudes. 

In this paper we will focus on unstably stratified flows, for which the approach is necessarily restricted to the early 
evolution of the interface. Although some of the nontrivial dynamic effects mentioned above are associated to the 
highly nonlinear regime, it may be useful to know within a controlled approximation to what extent those or other 
effects show up already at the early stages of nonlinear mode coupling. 

For the stably stratified case the weakly nonlinear analysis is obviously valid for long times since all mode amplitudes 
decay with time. Although this configuration may seem trivial, this is not the case in some situations, for instance 
when some external source systematically drives the interface out of its equilibrium state (planar or circular) . An 
example of this is the presence of noise sources, such as the quenched noise associated to a porous medium pH]. In 
a study of the long-wavelength, low-frequency scaling properties of the interface fluctuations, the knowledge of the 
lowest order nonlinear terms and their dependence on parameters such as viscosity contrast is crucial. In this context 
the weakly nonlinear expansion is the starting point of any Renormalization Group analysis of the relevant terms and 
of the fixed points of the problem. This line of research is clearly beyond the scope of this paper and will not be 
pursued here. 

The basic ideas of weakly nonlinear analysis developed here were first applied to viscous fingering by Miranda and 
Widom ]^9|,^0) to both channel and circular geometry (with fluid injection). The present work is in part an extension 
of those previous contributions in several directions, and in part a detailed study of selected particular situations to 
assess the validity and limitations of the approach. 

First of all, we provide a fully systematic methodology which may be carried out to arbitrary order. We explicitly 
calculate the results up to third order to include important situations for which the second order couplings, first 
discussed in Refs. ]^,^0), vanish identically, such as for A = or for configurations with up-down symmetry (which 
includes the physically relevant time dependent single finger solution with width one half). We give the explicit 
predictions both in real and Fourier space. 

In our formulation we also address the case of centrifugal forcing of Hele-Shaw flows |3l| . The experimental study 
of rotating Hele-Shaw flows has revealed a rich variety of new phenomena |27|,|3^,|33| . From a theoretical point of view, 
new classes of exact solutions with B = have been found [p4p5| | . The role of rotation in the possible suppression of 
finite time cusp singularities in the absence of surface tension has been discussed in |36| . From an experimental point 
of view, important differences in pattern morphology and new dynamical effects have been found for low viscosity 
contrasts |37j . It seems thus important to have this case included in the weakly nonlinear formalism. 

In addition, our study extends the earlier results of Miranda and Widom |2^,||(| with the discussion of the con- 
vergence of the weakly nonlinear analysis. We find the explicit exact criterion to assure uniform convergence of the 
series. Beyond that condition the series is asymptotic and different resummation schemes are also explored. The 
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different orders of approximation, including possible resummations, are carefully compared to exact solutions for the 
case of a single finger configuration. We find that in some cases the agreement even at relatively low orders is quite 
remarkable. 

The layout of the rest of the paper is as follows: in Section || we introduce the formalism. Section III deals with 



the derivation of the weakly nonlinear equations and their application to Hele-Shaw flows in channel geometry. The 
analysis of the radial geometry is carried out in Section [V , and Section |v| presents a numerical analysis of exact and 
approximate solutions. The main results and the conclusions are summarized in Section [v| . 

II. VORTEX SHEET FORMALISM 
A. Channel Geometry 

Let us first consider the Hele-Shaw problem in the channel geometry. We consider fluid 1 (viscosity fi±, density pi) 
below fluid 2 (p2, P2) (Fig. The z axis is perpendicular to the cell. A velocity Voo is imposed at infinity in the 
y direction. Gravity points from fluid 2 to fluid 1. The width of the cell is L, the gap between plates is b, and the 
surface tension between the fluids is a. 

The equations of motion for the interface and the boundary conditions are well known Q. Here we will use the 
formulation of Tryggvason and Aref We introduce the velocity U = (ui + tl%)/2 as the mean of the two limiting 
values of the velocities from both sides of the interface (iti, it 2) at a given point. This velocity U can be expressed in 
terms of the vortex sheet distribution 7 at the interface as: 

where: 

7 = 2A(U-s) + 2C , y-s + 2L>K s , (2) 
with s the arclength, k the curvature, and 

A = ; 7 C =Tv i r+AVoo, D — —— ■ -, 7= {u 1 -u 2 )-s. 3 

P2+P1 12(^2+ Mi) 12(M2 + Mi) 

For the purposes of this work it is convenient to rewrite these equations in terms of the interface height h(x), following 
p8[. The equations read: 



u( X t) ±p r w)-k*),*-*) , (x , w (4) 



7 = 2Dk x + 2Ch x + 2AU ■ (1, h x ), (5) 

where: 



7 = 7^1+^1, K " (1 + ^)3/2 - (6) 

The dependence of h and 7 on time is not written explicitly 

To complete the definition of the moving boundary problem the continuity of the normal velocity at the interface 
is required. This means that the velocity in the y direction, dh/dt, projected along the normal direction, is equal to 
the normal component of the average velocity of the interface, U ■ n: 

^=Uy- U,h x . (7) 



In Section HI we will consider the interface in the comoving frame (h — 0) and will look for the equation of evolution 
of h[x, t). 
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B. Radial geometry 



We proceed to write the corresponding equations for a circular cell rotating with angular velocity Q. The initial 
condition is R = Rq (constant) and we label the outer (inner) fluid as the fluid 1 (2) which have known viscosities pi, 
p 2 and densities p\, p 2 (Fig. ||). 

If we perform a change to polar coordinates in ([j]) , and from arclength s to angle <fi, we get for the mean velocity: 

u, = u(m -f=±p r ;^-^> (8) 

27r JO r l+ r 2~ 2r l r 2 COS(0 2 - Pi) 
1 f * T\T2 — r\ cos( 



* ' 2tt J rf + r|-2r 1 r 2 co 8 (0 2 -^ 1 ) RV2J ^' v; 



where 7 = y/l + (r<p/r) 2 (ui — w 2 ) ■ s, and we have used the notation r(<pi,t) = r%, r(0 2 ,i) — J" 2 . 

In the presence of sinks or sources the velocities u\ and w 2 which define U include only the solenoidal part of the 
total velocity field. For this reason, in the presence of injection (Q > 0) or withdrawal (Q < 0) of the inner fluid, Eqs. 
(||) and (g) must be supplemented with the corresponding irrotational part of the velocity field. 

In order to obtain the expression for the vorticity as a function of Uj, we use the local equations known for this 
problem p7pl: 



v pi = - 

and the boundary conditions 



V Pt = I Ui + ) + n 2 Pl rf, i = 1, 2 (10) 



(!) ( 2 ) Mil 



After some algebra, we can write an expression for the vorticity as a function of U: 



b 2 2a -> . . „ Q 6 2 2fi 2 (p 2 -Pi) 

12r(pi+p 2 ) 2nr z 12 pi + p 2 



with 

(r 2 + 2r 2 - rr u ) 



A= P2_P1 (13) 



(r 2 + r|) 2 P1+P2 

The last step is to derive the equation of the continuity of the normal velocity. Now the projection of the radial 
velocity along the normal direction n has two contributions: the solenoidal part of the average velocity, U, projected 
along h, and the irrotational part of the mean interface velocity, Ui rro t'. 

po=0 n + 0,^-n - l^-a^+A. (14) 

This completes the system of equations for the two geometries considered. For an interface with no overhangs, Eqs. 
(I), and (Q) are the starting point for the generalized weakly nonlinear analysis in the rectangular cell. Equations 
(3), @, (fl2"l), and ( fbil) play the same role for the analysis in the rotating circular cell. 



III. SYSTEMATIC WEAKLY NONLINEAR ANALYSIS. CHANNEL GEOMETRY 

Our goal in this section is to introduce a systematic method to derive an evolution equation of the interface in real 
space to a given order in nonlinear couplings, in the channel geometry. The different orders of mode couplings will 
be ordered as powers of a "book-keeping" parameter e, to be defined below. The evolution of the interface will thus 
take the form: 

^ = F[h] + eG[h] + e 2 I[h] + ..., (15) 
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where F[h], G[h], . . . are nonlocal operators on the function h(x, t), including nonlinearities of order n + 1 in the term 
of order e™. The small parameter e is denned as the ratio of two lengths, e = w/L. We take w as a measure of 
the characteristic scale of variation of the interface h(x), while L is either the width of the cell or, alternatively, the 
characteristic scale of variation in the x direction. The weakly nonlinear regime is defined by the condition w -C L. 

The order e° in Eq. ( [l5| ) corresponds to the linearized equation. The order e 1 , when written in Fourier space, 
corresponds to the result of Miranda and Widom |2{J. Here we will perform the explicit calculation to one order 
higher in e (up to I\h]) which is the leading nonlinear contribution in several important cases, such as the one 
discussed in Section HI D . The systematics is presented in the channel geometry for simplicity. 



A. Dimensionless equations, expansion and convergence 

We scale the interface height with w, the coordinate x with L, and the time with L/C, where the velocity C has 
been defined in (H). 

The equations (jl), (§), ©, and (Q) become: 



*f(M) = -p r {£[ v- hix)hx - x,) 2 ^(z<) d *i, 

2n J i r lui^I 



(x - x') 2 ll + e 2 



h(x') — h(x) 
x' —x 



(16) 



7 = 2Bk x + 2sh x + 2AU-(l,eh x ), (17) 

where 

ab 2 sh xx 
= 12CL 2 (^i+/x 2 )' K " (l + e 2 ^)3/2' ( 18 ) 

and 

dh 1 



dt e 



Uy - U & h x . (19) 



The starting point of our approach is an expansion of U, 7 and n in powers of e. Notice that the term 
\h(x') — h(x)] 2 [x' — x) 2 between curly brackets in ( |l6| ) is bounded provided that h x does not diverge. Equation 
( |l9| ) thus takes the form: 

f = ~ £ U^ + - h xU f + eiuf - h x U^) +..., (20) 

Before pursuing the calculation in detail, let us point out that the e-expansion in Eq.(pO|) has a finite radius of 
convergence. This is guaranteed by the properties of uniform convergence of both the expansion of the denominator 



in Eq.(16) and the curvature. These properties allow us to commute the expansion with the integral in Eq.(|i6[) and 
yield a convergent series of the form ( p0| ) . The radius of convergence of the expansion of the denominator of Eq. ( |l6| ) 
is given by the condition e 2 [h(x') — h(x)] 2 [x' — x} 2 < 1, while the convergence of the curvature expansion is assured 
by the condition e 2 h 2 < 1. In the nonscaled original variables, the two conditions for convergence coincide, and read 
MAX(\h x \) < 1. If \h x \ < 1 in the whole domain of integration, then the e-expansion converges. If the condition 
is not fulfilled in some intervals, then Eq.(|20|) is an asymptotic expansion. Even in this case, the expansion contains 
useful information about the original problem. 

An interesting case is A = which makes the vorticity independent of U in Eq.([l6|), so that Eq.(|l9|) becomes a 
closed equation for h(x, t). Then, in Eq.(po|) it is easy to show that = when n is even and = when n is 
odd, a property that makes the even power terms of the expansion in Eq.(po|) to vanish. 

B. First and second order expansion (e° and e 1 ) 

Following the scheme introduced in the previous section we obtain from Eq.(|l^) that: 
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» 2ir J_ x x-x' y 2ir x-x' K ' 



Since -y(°) = 2Bk x °^ + 2AU^ — we get f/J — U- ^ — 0. With this result we can study the first order term in the 
vorticity equation: 



7 «(z) = 25k« + 2h x + 2A\jf + 2Ah x uf ] . (23) 
Taking into account the definition of the Hilbert Transform: 

H[f{x>)] = \ P £j ^ d *'> ( 24 ) 



Eq.(p0|) up to order e° reads: 

dh 



-H[Bh x , x , x> +h x ,]. (25) 



The linear operator F[h] in Eq.(|T^) thus reads explicitly: 



F[h] = -P [ + °° ' n, '°; '■ <!,'. (26) 



Writing h(x,t) as a superposition of Fourier modes in Eq.(|25|) we recover the linear dispersion relation: 

S k (t) 



s k (t) 



X(k) =| k I (1 - Br). (27) 



We will take k as an integer but we should keep in mind that, upon restoring dimensions, k should become (2-7T / L)n, 
with n integer. 

Let us pursue the systematics of the method by computing the next order in e in Eq. ([l^) : 

^ = U^+sU^=F[h}+eG[h]. (28) 

(2) 

The computation of U$ requires an expression of the vorticity up to second order. This includes the evaluation of 

(2) 

U~ , which must be computed from the first order term of the vorticity. We get: 

Uf = H [{h{x')f{x')) x ,\ - h(x)H[f x ,(x% (29) 
with f{x) = (Bh xx + h) x = 7 (1) (aOA and: 

Uf ] =A{H [h x ,H [/(*")]] + H [h(x')H[f x „ (x")}] + {hf) x } . (30) 
Taking the Fourier transform, we obtain: 

+00 

S k = X(k)5 k +sA\k\ [l-sgn(ks)]X(s)5 s (t)S k - s (t), (31) 

s— — 00 

which coincides with the result of Miranda and Widom in Rcf. l29| . 
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C. Third order expansion (e 2 ) 



The expansion to order e 2 is necessary to account for the lowest order nonlincarities in the case of zero vi scosit y 
contrast, and for other relevant situations such as the time dependent Saffman-Taylor finger solutions (section [II D| ). 
We now have: 



f = O(e ) + O(e 1 ) + £ 2 (C/f -Wtf>) + - 



We have already computed in Eq. (|29|). On the other hand 

2 -l; v~ /J ' 2 7T P 



(32) 



(33) 



Integrating twice by parts, the last integral can also be written as a Hilbert Transform. After some algebra we obtain 
the explicit form of the operator I[h] containing the cubic nonlinearities in Eq.([Lq), which reads: 



with: 



and: 



I[h] = + \H [(f(x')[h(x) - h{x')f) ] 

~h x H [(f(x')[h(x') - h{x)]) x ,\ + V[h, A] 



V[h,A] = A 2 H[h(x')H{T x „{x")\ + h x ,{x')H{T(x")\\+A 2 (hT) x , 



-TTW^h rr(l) _ 7 (2) (z) 



l(x) = B (h xx h 2 x ) , t(x) = U x z) +h x U\ 



2A 



The same result in Fourier space takes the form: 

+00 



8 k {t) = 0[l,e) + e 2 J2 W-ih- 



A 2 T{k, s, I) - -BY(k, s, I) + W(k, s, I) 



(34) 



(35) 



(36) 



(37) 



with: 



T(k, s, I) =| k 1 1 a I A(l) [1 - sgn(fc S )] [1 - sgn(J S )] , 



(38) 



F(fc,s,Z) =| k I / 2 (s - 0(fc - s), 



(39) 



W(fc,s,0 



Z \ fc 2 

Z ( - + k — s J — sk sgn(Zs) + — sgn(fcZ) 



A(i). 



(40) 



It is clear that the viscosity contrast in Eq. (|37j) appears squared because of the reflection symmetry (the simultaneous 
change A — > —A and ft, — > — /i is a dynamical symmetry of the problem). Symmetry reasons alone, however, do not 
allow to discard a three mode coupling contribution when A = 0. We see from our calculation that a three mode 
coupling is indeed present independently of A. 

Following this scheme, the fourth order will carry a contribution proportional to A, and another proportional to 
A 3 , for symmetry reasons. The fifth order will carry a contribution independent of A and two others, proportional to 
A 2 and A 4 , respectively. This scheme will continue for subsequent orders. 



D. Analysis of the time dependent single finger solution 

In this section we perform a detailed analysis of the weakly nonlinear expansion in cases where exact solutions are 
known, namely single finger configurations with B = 0. This allows to study how exact properties of solutions show 
up at the different orders, particularly concerning the role of viscosity contrast A. 
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At this point it is worth recalling that the case A — 1 allows for a continuum of Saffman-Taylor finger solutions 
corresponding to different finger widths. A continuum of time dependent exact solutions, leading to those stationary 
states, is also known for A = 1. However, for A ^ 1 only the one of width A = | remains a solution, A b eing the ratio 
of the finger width to the width of the channel. This result, which has been recently addressed in Ref. [g6| although 
it was first discovered in Ref. p5|] , points out an intriguing connection between the width selection problem and the 
dynamical role of viscosity contrast. Here we will analyze the interplay between A and A in the early nonlinear regime, 
and elucidate at what stage of the nonlinear dynamics does the viscosity contrast 4 / 1 prevents the possibility of 
having A / 1/2. 

From now on we consider L = 2tt, B = 0, C = 1, and we take as scaling velocity (1 — A)C = (1 — A) = Conformal 
mapping techniques enable us to write the single finger solution of this problem in the form ]25|,^6| : 

f(w,t) = -Inw + d(t) + ?/ln(l - a(t)w), (41) 

where f(w, t) = y + ix, is an analytic function inside the unit disk in the w-complex plane, which maps that disk 
into the physical region occupied by the more viscous fluid. The interface is obtained in a parametric form by setting 
w = e lB . The functions a(t), d(t) verify: 

4*) = ^fS-?f). (42) 



(43) 



2-77 V 2 a(t) 

a(t) 



«(*) 2 + rKry-2) +?? (2-77)i±g$ 



To obtain an expression for h(x) in the weakly nonlinear regime of the evolution, a(t) and d(t) are expanded in powers 
of a small parameter v. 

d{t) = d {a) {t) + vd {1) (t) + v 2 d (2 \t) + --- , a(t) = va (0 \t) + v 2 a {1) {t). (44) 
Introducing these expansions in Eqs.(f42|) and (|4^) we obtain: 

$o) =l p-)=0, d< 2 > = ^(a<°>) a (t), d< 3 > =r ? 3 a(°)a«, ••• (45) 

d< > = V°>, dt^V, d( 2 ) = ^( a ( 2 )- 77 (2- 77 )(a(°)) 3 ), ■•• (46) 
which will be useful later. From ([ll]) and ( [l2| ) we obtain: 

y = h= -urja^ cos 6 -v 2 (tja^ cos 6+ ^(a^fcos2d-d^) , , 

-v* (rja^ cos 9 + -qa^a^ cos 29 + §(a(°)) 3 cos 30 - d^) + 0(j/ 4 ), ( ' 



x = — 9 — vqa^ sin( 



t/ 2 (^ (1) sin0+|(a (o) ) 2 sin20) + ©(t/ 3 ). (48) 

This last equation can be inverted in a systematic way to get: 

9 = -x + viia {0) sinx + v 2 (r/a^ smx + |(1 - 7?)(a (0) ) 2 sin 2^ + 0(v 3 ). (49) 

The relation between 9 and a; can be inverted in Eq. (f47|) and the cosine functions expanded. We obtain in this way 
the following expression for h(x, t): 

h(x, t) = —vrjoc-®* cos x — v 2 r\ (a^ cos x + (a^ ) 2 cos 2x) + 

z/ 3 7y {(77 - cos2x +(f{v- 2)(a(°)) 3 - a { ^) cosx+ (50) 

(f (2 -77) - i) (a(°)) 3 cos3.T} +C(^ 4 ). 

In order to follow the scheme developed in the previous section, we must measure h(x, t) in units of its characteristic 
amplitude v. In this way, the previous equation (E0) takes the form: 
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h = h<® +// 2 /i (2) +0(v 3 ), 



(51) 



where h(x,t) represents now h{_x,t)/v, and the small parameter v is directly comparable to the small parameter e 
of the previous section. The expression for h(x,t) can be regarded as a series of modes of decreasing amplitude in 
our expansion ([l5|) (written in the units of this problem), and matching the corresponding powers in either e or v we 
obtain: 

m 

to order e°. 

This identity can be verified (independently of A) by substituting the expression of hW and using Eqs.(|4"6|) and 
(Eq) for the left and right hand side respectively. At order e we have: 



dt 



= Vjr[ft(i)] + G[M°),ftW]). (53) 



The second term of the right hand side gives no contribution, and we obtain an interesting result: the equation at 
order e is verified independently of r\ and A. Hence, we cannot establish the difference between A — 1 (compatible 
with any 77) and A 7^ 1 (compatible only with r\ = 1) at this order. 
The difference between these two situations arises at order e 2 : 

^ = | {F[h^\ + G[h<®,h<»] + L[A<°U<°>, A<°>]) . (54) 

Since the left hand side of the equation involves only the modes cos a;, cos 2x, and cos3x, the right hand side of the 
equation must include these same modes with the same coefficients. The coefficients for cos2x and cos3x are easily 
matched. For cos a: the left hand side reads: 

^^( a (°)) 3 -^ 2 ))cos,, (55) 
where we have used Eq. (|46j), and the right hand side reads: 
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ll a (°)f _ IL a (A cosx + A 71 £— i£ ( a (0))3 cosx + ^{a^f , 
2 J 4 8 



(56) 



Hence, the requirement for matching the coefficients on both sides is: 



^(l-»7)=^(l-fj) > (57) 

which is always true if 77 = 1 (A = 1/2), but for ry 7^ 1 (A 7^ 1/2) it requires A = 1. This result shows that the nontrivial 
relationship between A and A which is known from exact solutions is already manifest at the early nonlinear stages of 
the dynamics. This clearly illustrates the potential usefulness of the weakly nonlinear expansion at a purely analytical 
level, in that a dynamical property of the problem which must be satisfied at all orders (in this case the incompatibility 
of A 7^ 1 and A 7^ 1/2) may be detected at a low order in the expansion, without having to know the exact solution 
to all orders. 

If we pursue the expansion to higher orders, the general expression for hM 1 ' takes the form: 

n+i 

=J2Pi n) ( t ) C0S ( kx ): (58) 
k=l 

where each coefficient is a function of a(t). The even modes of the solution with 77 = 1 have zero coefficients because 
they must remain invariant under a change of sign and translation by n. 

So far in this section we have restricted ourselves to B = in order to compare with exact results. However, we 
can carry out the analysis also for B^O. It can be shown that the structure of the expression (|5^) is preserved in 

this case. The coefficients (3^ (t) are obtained as solutions of a set of differential equations which contain the surface 
tension B. For instance, to third order we have: 







h(x,t) = (/^ 0) +sp[ 1} +£ 2 M 2) ) cosx+ Ltfp +£ 2 /?2 2) ) cos2x + £ 2 ^ 2) cos3x, (59) 

which, introduced in Eq.(|l5"|), leads to the following equations: 

pf ] = V (l - B)(3f\ /3« - F (l - = 2F (1 - 45)^ 1) , 

= 2%(1 - 4S)M 2) , /f } = 3%(1 " 95)/f } - ^(/3< 0) ) 3 , (60) 
/3f ' = V (l - B)0™ + AV (1 - B)^^ - y (2 - Si?)^) 3 , 

once the dimensions are reintroduced. In this way we can also see how surface tension perturbs the dynamics at 
the weakly nonlinear regime. Clearly, at these early stages of the nonlinear evolution there is no sign of the singular 
perturbation character of surface tension, which, at the late stages of the evolution will be responsible for the selection 
of the steady state Q . 



IV. APPLICATION OF THE METHOD TO THE ROTATING HELE-SHAW CELL 



A. Dimensionless equations in the radial geometry 



The e expansion in the channel geometry originated in the different scaling of lengths in the x and y directions. 
To perform an equivalent expansion in the radial geometry we first split the function r(</), t) in two parts, namely the 
radius of the unperturbed circle and the departure from this circle: 



K&i) \/i? 2 + ^ + r(<M) = R(t)+r(<j>,t). ((ill 



7T 



The largest length scale in the problem is given by R(t), which is only a constant in the absence of injection. For 
finite injection rate, R(t) defines an evolving (unstable) solution. The proper scaling of r((j>,t), which will define the 
weakly nonlinear dynamics now, is not the naive extension of the scaling defined in the channel geometry (namely 
£ = w/R, w being the typical departure from circularity). In fact, this scaling would mix different orders of the 
weakly nonlinear analysis. The reason for this was already pointed out by Miranda and Widom |30[ |: in the radial 
geometry, mass conservation implies that the zero mode (which in the channel geometry is decoupled from the rest 
and drops out of the formulation) has a higher order nonzero amplitude, since it must satisfy a mass conservation 
relation which, to lowest order, reads: 

Following this observation, we split the deviation from R in two terms, so that the radial function reads r = R+f+ro, 
where r$ represents the zero mode (an expansion or contraction of the circle), and f the other modes. We scale them 
in the form f — » wf, ro — » (w 2 / R)ro, and write the position of the interface r{<f>) as: 

r((j>) -> R(l + er(4>)), (63) 

where: 

r(4>) = r{4>) +er . (64) 

To simplify the notation, we have dropped out the time dependence. 
The velocities will be scaled by: 

v ° - t4t" i^rir 1 - t^ 2 - <*>*) • ( 65 ) 

implying that time will be scaled by R/Vq. Once equations (||), @, ([l2|), and ( |4j ) are made dimensionless, we have: 

7 = 2B-^- + 2AU ■ 1) + 2e (j-^ d) r>, (66) 

1 + er 1 + er \{l + er) z J 
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for the vorticity, where: 



B 



12 ( m + fi 2 )V R 21 



C 



QA 
2irRV ' 



D 



b 2 il 2 ApR 
12 +^2 )Vb' 



(1 + er) 2 + 2e 2 r 2 - e (l + er)r^ 



(1 + er) 2 + £ 2 r 2 



- A*i 

M2 + Ml 1 



A/? = p 2 - pi- 



le?) 



For the integrals: 



1 + 2er 2 + e 2 r| 



o tan 



92—91 



[1 + e/(^i, ^ a )+e 2 ^i,^ 2 )] 



-7(^2, *)#2 



(68) 



2tt 



1 + e/(0i, 02 ) + £ 2 rir 2 - (1 + 2er 2 + e 2 rf ) cos(0 2 - 0i) _ 



2 shT 



2 / <fe-<; 



^) [l + £/(0 1 ,02)+£ 2 5 (0l,0 2 )] 



7(0 2 )d02, 



(69) 



with: 



/(01,02) = Tl +r 2 , 5(01,02) 



and for the evolution of the deviation from i?(£): 



(ri - r 2 ) 2 



4sin' 



/ 92-91 \ 

V 2 y 



rif2, 



(70) 



£ ^ £ eft '' £ l + sr 2nRV Q 



1 + er 



(1 - e 2 r ) 



(71) 



Eq. ([7l|) above is the only one which is sensitive to the different possible scalings of the deviation from circularity, 
due to the time derivatives. The choice we propose is the only one which is truly systematic. Other possibilities are 
not incorrect but will mix different orders of the weakly nonlinear expansion at a given order in e. 



B. The linear dispersion relation 

Our goal now is to obtain the linear dispersion relation and the leading weakly nonlinear corrections, in both real 
and Fourier space, and to discuss the interplay between rotation and injection at the early stages of the instability. 
We will use the following definitions for the average of a function and for the Hilbert transform in the unit circle: 

1= hfJ fm ^ H[m] = h p C f ^ coi (^) 64,1 ■ (72) 

The zero order of the vorticity and the two components of the velocity are: 

7 (0) = 2AC/f, ^f = ?> ^ O) = ^[7 (O3 (0')], (73) 

and thus: 

7<°> = A^. (74) 

This is different from its counterpart in the channel geometry. Here the equations for the vorticity at consecutive 
orders require the knowledge of the average vorticity. The equations are solved by averaging on both sides. For 
example, at zero order, yl / 1 will lead to = 0, but for A = 1, 7W can be any constant. The presence of 
an arbitrary constant is not a problem, since it is eliminated at each order in the equation for the evolution of the 
interface. 

Expanding the equations up to first order, and using 7^ = 0, the linearized equation for the evolution of the 
interface deviation reads: 
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with: 



7« = -2B(f + r^) 4 , + 2{C-D)P 4> . (76) 
We can now perform a Fourier transform of the equations and, after reintroducing the adequate dimensions, we find: 

(b 2 » 2 A P \ 6 2 , Q 

This is the linear dispersion relation found in pjj] . 

We now return to the question of the different scaling alternatives introduced at the end of Section [V A . Had we 
used the scaling e = w/R 2 (t), the term Q/(2ttR 2 ) in ( |77| ) would have disappeared before reintroducing the dimensions. 
On the contrary, for e = w this term becomes multiplied by a factor 2. These changes in the dimensionless linear 
dispersion relation come from the fact that a mode S n in the scaling adopted in (^) and (Q) becomes S n /R in the 
first of the scaling alternatives above. 

The term Q/(2ttR 2 ) at the end of Eq. ( |77j ) is stabilizing for positive injection rate (Q > 0) and destabilizing for 
Q < 0. In a sense this is a purely geometric effect, which contributes to the growth or decay of modes due to the 
expansion or contraction of the base state. As pointed out in J3(| the term Q / (2itR 2 ), which can be scaled out, must 
be distinguished from the rest, which do describe the intrinsic instability of the problem. Disregarding that last term, 
it is then clear that injection and rotation (CVq and DVq respectively) play an equivalent role in the linear dispersion 
relation, since both appear multiplying n. By choosing the proper viscosity and density contrasts, they can produce 
exactly the same stabilizing or destabilizing effect. This is the counterpart, in the radial geometry, of the equivalence 
between the roles of injection and gravity in the channel geometry. In the channel geometry, however, the equivalence 
is exact and can thus be extended to all the nonlinear evolution (in the appropriate reference frame). This is not the 
case in radial geometry, as shown in the next section. 

C. Weakly nonlinear theory with rotation 

The purpose of this section is to derive the leading order nonlinear contributions for the general problem with 
both rotation and injection. We first recall that mass conservation related the zero mode to the other modes. In the 
original nonscaled variables this reads: 

^ = -2^°-^5> fc|2A(fc) - ( 78 ) 

To reproduce this result and the nonlinear couplings for the other modes we start with the equation for the interface 
at order e 2 : 

* + ^.^ + .(^ + _|_^), (79) 

where we have used that = E/i = U- ^ = 7^ = 7^°) = 0. The velocity Up can be obtained from the expansion 
oiUf. 

Uf ] = i#[7 (2) (0')] + - m)l (1 H<P% (80) 



where: 



7 (2) = 2Bk {2) - 2Sf K (1) + 2Af 4> uf ) + 2AU (2) - 4Cff> (81) 



K ( 2 ) and U\ > can be computed from (|67]) and ([39]) respectively. The latter involves only 7W. We obtain finally: 

J + = O(e ) + e (h [Bnf - ^(r%] + fff [ V ] + AH[H[{rs r ) r ]] 

- AH [f(#)H[a rr ] + rpH[sp,]] + ^r 2 ) , 
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with: 



k (2) = r 2 + -| + 2ff^, S0 = (B(f + + (D - C)f) q 



(83) 



By introducing a superposition of Fourier modes in f we directly recover Eq. (|78|) . For the other modes (k ^ 0) we 
get the result: 



5 n = \(n)5 n + S k S n -k{F(k,n)+S{k,n) + X(k)J(k,n)), 



where: 



F(k,n) 



R 



QA (\ 



2irR 2 V 2 



sgn(fcn) 



b 2 n 2 A P 



24 p x + p 2 



(84) 



(85) 



S(k, n) 



|nj 
i? 



12i? 3 /ii 



■(n + 3Jfe) 



(86) 



J(fc,n) 



i? 



[A | n | (1 - sgn(fcn)) + 1] 



(87) 



For the particular case O = 0, the expression above reproduces the result of Miranda and Widom 30 1. We find that 
the presence of rotation does not change the formal structure of the equations, but introduces new terms. 

To emphasize the roles of injection and rotation, we define the quantity H(k,n) as the part of the coupling matrix 
in the r.h.s. of Eq.(|84|) which contains Q and/or fl. This quantity reads: 



QA b 2 n 2 Ap 



2nR 2 12 /ii + [it 



J(fc, n) | k 



b 2 n 2 Ap QA 



12 fn + p 2 2nR 2 



2nR 3 ' 



We observe that experimental parameters occur only in two groups, namely: 



b 2 n 2 A P 

12^ 



pi 



Q 



QA 
2nR 2 ' 



(89) 



Both of them turn out to multiply the same functions of the wave numbers, except for changes in sign. However, 
notice that the relative sign of f2 and Q is different in the linear part and in the leading nonlinear part of the dynamics. 
This clearly shows that the effects of injection and rotation cannot be interchanged by simple changes in parameters, 
and that the intrinsic new dimensionless parameter related to rotation already shows up at the early nonlinear regime. 

It is also remarkable that, for Q ^ 0, R(t) is time dependent, and the relative weight of fl and Q changes with 
time. This may have remarkable consequences in the highly nonlinear regime. For instance, rotation can prevent the 
formation of finite time singularities in the case with zero surface tension J3(| . 

Finally, notice that at nonlinear order there is also a geometrical term, Qj (2irR 3 ), independent of n. This term 
cannot be scaled out by the same procedure that eliminated its counterpart at linear order. 



V. NUMERICAL STUDY OF AN EXACT SOLUTION AND ITS WEAKLY NONLINEAR 

APPROXIMATION 

The purpose of this section is to put the weakly nonlinear analysis to the test, as a quantitative approximation. We 
will check it against an exact time dependent solution of the case B = which evolves towards the Saffman-Taylor 
finger of width A = 1/2 in the channel geometry. This will allow us to check how fast is the convergence of the 
expansion and how accurate it may be even beyond its radius of convergence, when it becomes an asymptotic series. 

We define F as the ratio between the maximum height of the interface (at the finger tip) and half the width of the 
cell (in our case L = 2tt) . F is a dimensionless amplitude which is proportional to e and thus measures how deep the 
system is into the nonlinear regime. Furhermore, since the system is described by the evolution of a curve, it may 
also be convenient to have a more global characterization of its configuration in order to compare the exact result and 
the different approximations. We propose the use of the flux $, defined as the total amount of fluid 1 per unit length 
and unit time flowing across the horizontal line located at the mean interface position fig]). 
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A. The time dependent exact solution with B = and A = 1/2 



An explicit solution of the problem without surface tension (B = 0) and valid for any viscosity contrast A, describing 
the growth of a single finger which asymptotically fills half of the channel (A = 1/2), can be obtained from Eqs.(^Tj), 
@, and (p|). This reads: 



f(w,t) =-lnw + V o t + a(0)+]n y/l + (b(0) 2 - 1) exp(2V t) + 



w 



with a(0) and 6(0) as initial constants, and 6(0) 3> 1. Since this solution is valid for any A, particularly A 
terms of the weakly nonlinear equation (15) depending on A will necessarily have no contribution. 
We will take the following values for the parameters: 



v = ~, 6(0) = ~ = iooo, 



<z(0) = hxv + — . 



(90) 



0, all 



(91) 



Fig. U shows the evolution of the interface for this solution. When F ~ 1 the finger shape is indistinguishable from 
the stationary Saffman-Taylor solution (Fig. ^), except near the points x — ir/2 and x — 3tt/2 when the asymptotes 
of the finger must develop at infinite time. When F = 1/tt, the slope of the in terfac e becomes 1 at these two points 
and the series (|l5| ) loses its convergence, according to the discussion of Section III A . 

In Fig. U we present the amplitude of the Fourier modes of the exact solution. It is clear that relatively few modes 
suffice to accurately reproduce the stationary finger shape in the tip region. 



B. Expansion of the exact solution 



In Section [IIP we expanded the exact solution (|4lJ) in a small parameter v. This expansion revealed a hierarchy 
of amplitudes of the different modes for an initial condition close to the planar interface. The solution ([)0]) can be 
expanded similarly and yields: 



In agreement with (|50|), and using the same parameters as in the previous section, we obtain up to third order: 



h(x, t) = ve 1 / 2 cos x + v 



L*/a _ \e 3t l 2 \ cosx- i e 3t / 2 cos3z 
2 8/ 24 



(92) 



(93) 



The exact flux $ can be easily computed from the stream function, which in its turn can be related to the form of 
the mapping Eq.(p0|) (see for instance Ref. In Fig. || we compare the flux $ of the approximate result above 

with the exact result coming from (0). We clearly observe that the first correction to the linear regime improves 
the flux (to a 5 % accuracy) from F ~ 0.12 to F ~ 0.17. The amplitude of the mode k = 1 (Fig. ||) is fairly well 
reproduced by the linear approximation up to F ~ 0.20, and the correction of order is 3 gives a value for the mode 
k = 3 reasonably good up to F ~ 0.20 (Fig. |^) and improves the mode k = 1 almost until F ~ 0.25. Above these 
values of F the deviation from the exact interface is exponential. 

Notice that the expansion up to v 3 considered here does not account for the complete three mode coupling of 
Eq.(|37|) which contains part of the higher orders in v. This explains why going from order v (linear) to order v 3 
improves only slightly the shape of the finger, the flux <E>, and the amplitude of the mode k = 1. 



C. Weakly nonlinear expansion 

In this section we study the mode coupling equation (^) . We solve this equation in cases where only modes k = 1 
and k = 1,3 are present, and compare the result with both the exact solution and the approximation to order v 3 
introduced in the two previous sections. 

We consider the initial condition /?i(0) = e = 0.001 for the first mode and zero for the rest of modes. We use Eq. 
( |37| ) recalling that, since /3j represent the amplitudes of the cosine functions, Si — Pi/2. The result is that the modes 
k^ 1 remain zero and: 
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= £ exp(0.5t) m 

v/l + 0.25e 2 (exp(i)-l)' 

where we see that the mode k = 1 saturates to a finite value as t — > oo. Returning to Fig. || and 0, we see that the 
flux $ approximates the exact value to a 5 % uncertainty) until F ~ 0.23. Fig. || shows that the amplitude of the 
first mode starts deviating significantly from the exact solution around F ~ 0.40. 

Next, we add the third mode to the equation ( p7| ) with the initial condition /3j = e + 3e 3 /8 and /?3 = e 3 /24. The 
set of equations to be solved numerically is: 

fa = \fa - \pl - + (95) 

3 3 27 

03 = ~ -fifo - -jPl (96) 

As shown in Figs. ^| and || the flux does not improve significantly, while the amplitude of the mode k = 3 remains 
acceptable almost until F ~ 0.25, and improves the previous result by more than a 50 % almost until F ~ 0.33. 



Beyond this value of F the third mode falls to zero due to the second term in the r.h.s. of Eq.(96), proportional to 
f3fp3 but negative. This is a signature that higher order terms, for instance of order f3f03, become important when 
[3i is of order one. 

From this study we can conclude that using only two modes and the lowest order nonlinear correction for this case 
(involving three mode couplings), both the shape of the interface and the total flux are reasonably well described 
within the whole range of convergence of the series (up to F ~ 1/tt). This indicates that the weakly nonlinear 
description of the problem works remarkably well at the quantitative level, at least in some physically relevant cases. 
To what extent this conclusion holds for more complicated situations, for instance those involving competition of 
different fingers, remains an open question. 



D. Partial resummation 



According to Eq. (£37|), we can write a closed equation for the mode k = 1 which is correct up to third order 
couplings, of the form: 

fa = l/h-~ftfa, (97) 

where the mode k = 3 has been set initially to zero. 

The presence of 0i in the right hand side gives rise to terms of order higher than j3f. Thus, solving this equation 
amounts to doing a partial resummation of these higher order terms. This kind of closed differential equation was 
already obtained in Ref. p9|| . In our approach this differential equation is obtained in a systematic way by identifying 
the terms A(n)/3„ and replacing them for n . For example, in Eq. ( B7j ) the term T(k, s, I) could have been partially 
resummated in the two mode coupling differential equation. 

The solution of Eq. ( |97| ) yields the transcendental form: 

t = 2 ]n(|)+i03?- e a ). (98) 

A numerical computation of the flux based on this equation shows a dramatic improvement up to F ~ 1/tt. It is 
surprising and remarkable that the equation for a single mode, including the resummation suggested by the first 
nonlinear correction, reproduces the exact solution to a great accuracy in the full range of convergence and defines 
an excellent approximation further inside the nonlinear regime. This suggests that this kind of resummation is not 
arbitrary, and might have a deeper physical justification. 

In addition, we obtain an amplitude for the mode k = 1 within 10 % accuracy up to F ^ 0.8, and it is still quite 
good until F = 1. Although the exact spectrum shows the increasing importance of the mode k = 3 at these values 
of F, the agreement between the exact result and this approximation for a single mode k = 1 is quite remarkable. 

We recall that the series converges up to F ~ 1/tt. Beyond this value of F, increasing the order of mode coupling 
in Eq. ( p7| ) does not guarantee that the approximation improves. Accordingly, there is an optimal finite order of 
approximation for a given value of e, as usual in asymptotic series. In this sense, our results above seem to indicate 
that Eq.(|97j) is the optimal approximation for the mode k = 1 in the asymptotic (nonconvergent) region (i.e., for 
moderately long times). 
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In the same spirit of Eq. (|97|), we can also derive a closed system of differential equations which includes third 
order couplings for the modes k = 1, 3. Its solution shows no qualitative differences with the case of the single mode 
k = 1 above, as far as the flux and the amplitude of the first mode are concerned. The improvement in the amplitude 
of the third mode and consequently in the shape of the finger does not reach the high values of F that the first mode 
reaches (F ~ 0.8), as we can see in Fig. |9[ In a short range of F the mode k = 3 is better than the case shown in the 
previous section, but deviates from the exact solution much earlier than the mode k = 1 for the same reasons than in 
the previous section. Specifically the mode k = 3 is extremely good until F ~ 0.25, starts deviating by more than a 
10 % at F ~ 0.33 and is about 50 % of the value at F ~ 0.40. 

We conclude that the weakly nonlinear formalism, apart from its nominal range of validity for very small amplitudes 
where it converges very fast to the exact dynamics, can describe regimes beyond those expected in principle with 
reasonable accuracy. In the case of the growth of a single finger, for instance, we have explicitly seen that simply two 
modes together with an appropriate partial resummation, yields a remarkably good approximation in the full range 
of convergence of the expansion (F < 1/tt). 

VI. CONCLUSION AND PERSPECTIVES 

We have developed a systematic scheme to derive the successive orders of mode couplings in a weakly nonlinear 
regime, adapted to the study of interfacial dynamics in Hele-Shaw flows. This has been done with full generality, 
including both the channel geometry driven by gravity and pressure, and the radial geometry with arbitrary injection 
and centrifugal driving, which includes the case of a rotating Hele-Shaw cell. The method could also be applied to 
even more general (position dependent) drivings. The formulation in real space has enabled us to address the issue of 
convergence of the mode coupling expansion. We have found that the exact condition for convergence in the channel 
geometry is \h x \ < 1 in every point of the interface. 

We have carried out the explicit derivation of nonlinear couplings up to third order in the case of channel geometry, 
in order to obtain the leading nonlinear contributions in cases where the second order contribution vanishes. These 
include the case of zero viscosity contrast and the time dependent single finger solution of width A = 1/2. 

On the analytical side, we have also illustrated the usefulness of the weakly nonlinear analysis in elucidating the 
role of the different parameters on the dynamics, in some examples. In the case of single finger solutions growing in 
a channel, we have shown how an exact nontrivial property of the problem (the relationship between A and viscosity 
contrast A for zero surface tension) can be extracted from an analysis order by order, without really knowing the exact 
solution. In the case of rotating Hele-Shaw flows, we have discussed the asymmetry between the roles of injection and 
rotation at the lowest nonlinear order, as opposed to what happens in the channel geometry, where the equivalence 
between gravity and injection driving is exact. 

On the numerical side we have checked the predictions of the weakly nonlinear analysis, at different orders, against 
exact solutions. We have found that the convergence is quite fast for small amplitudes. Furthermore, in the case of a 
single finger growing in a channel we have explicitly seen that the analysis to low orders (up to three mode couplings) 
yields a good description of the interface dynamics even close to the radius of convergence. With an appropriate 
resummation scheme, the prediction is relatively good even beyond that point. 

Some of the most interesting applications of the systematic weakly nonlinear approach have not been developed 
here since they deserve a separate, in-depth analysis. One of them is the study of the scaling properties of fluctuations 
in stably stratified Hele-Shaw flows with external noise pS| ]. A more direct application of our formalism which also 
deserves a separate analysis is the derivation of amplitude equations for the region near the instability threshold 
through a center manifold reduction. This would made contact with what is most commonly referred to as "weakly 
nonlinear" approach in the literature of pattern formation. A detailed study of this point with a careful discussion of 
the nature of the bifurcation and its implications will be presented elsewhere p9| . 
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FIG. 5. Amplitude of the Fourier modes of the exact solution at consecutive values of F. 




FIG. 6. Flux $ as a function of the factor F for different approximations (see text for details) . 



20 




FIG. 7. Relative error of the flux $ as a function of the factor F for different approximations (see text for details). 




FIG. 8. The first mode amplitude as a function of the factor F for different approximations (see text for details). 
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FIG. 9. The third mode amplitude as a function of the factor F for different approximations (see text for details). 
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